##             Rplots--Widower_Effects_in_Voter_Participation
##             ==============================================

## Author: Will Hobbs

## Table of Contents
## =================
## 1 Main (9 plots)
##     1.1 Preamble
##     1.2 Read files
##     1.3 Plot
## 2 Wide
##     2.1 Preamble
##     2.2 Read files
##     2.3 Plot
## 3 Gender comparison
##     3.1 Preamble
##     3.2 Read files
##     3.3 Plot
## 4 Ages
##     4.1 Preamble
##     4.2 Read files
##     4.3 Same voting history plot
##     4.4 Voting history comparison plot
##     4.5 Age and turnout
##         4.5.1 Treatment and age
##         4.5.2 Control and age
## 5 Voting History Comparison
##     5.1 Preamble
##     5.2 Read files
##     5.3 Plot
## 6 Recovery Plot
##     6.1 Preamble
##     6.2 Read files
##     6.3 Plot
## 7 Aggregate numbers
##     7.1 Preamble
##     7.2 Read files
##     7.3 Aggregate


## 1 Main (9 plots)
## ~~~~~~~~~~~~~~~~~
##
## 1.1 Preamble
## =============
  rm(list = ls())

  option.name <- c("All Registered Voters")

      x.title <- "Weeks Since Spousal Death"
      y.title <- "Proportional Widower Effect"

      election <- c("ss9", "gp10", "gg10")

      election.name <- c("Special Statewide 2009", "Gubernatorial Primary 2010", "Gubernatorial General 2010")

## 1.2 Read files
## ===============
  file <- NULL
         readfile <- NULL

          for (i in 1:length(option.name)) {

         for (j in 1:length(election)) {
            file[(i - 1)*3 + j] <- paste(option.name[i], " ", election[j], sep = '') # filenames no ext.
        }
        }
           readfile <- paste("output/", file, ".csv", sep = '') # filenames + dir. & ext.
           file <- gsub(" ", ".", file)

           loess.slice.1 <- paste(file, ".loess.1", sep = "")
           loess.slice.2 <- paste(file, ".loess.2", sep = "")

            for (k in 1:length(readfile)){
             assign(file[k], read.csv(readfile[k], head = TRUE))
             assign(loess.slice.1[k], get(file[k])[get(file[k])$week<=0,])
             assign(loess.slice.2[k], get(file[k])[get(file[k])$week>=0,])
           }

## 1.3 Plot
## =========
  for (i in 1:length(option.name)){
    direct <- paste('plots/Cohort -- ', option.name[i], sep = '')

    # prepare legend variables
    filecases<-round(mean(get(file[(i - 1)*3 + 1])$cases.in.week.matched))
    filecases<-paste("# cases/week ~=", filecases[1])
    filecontrols<-round(mean(get(file[(i - 1)*3 + 1])$controltocase.ratio))
    filecontrols<-paste("controls:cases ~=", filecontrols[1])
    filecontrols<-paste(filecontrols[1], ":1", sep = '')
    match.ratio<-round(mean(get(file[(i - 1)*3 + 1])$match.ratio), 2)
    match.ratio <-paste("match rate ~=", match.ratio[1])

    # prepare plot variables

    color <- c("goldenrod", "forestgreen", "purple4")
    span.par <- list(c(.3, .3), c(.5, .4), c(.3, .3))
    x.axes <- list(c(-91, 4), c(-31, 59), c(-15, 80))
    y.axes <- c(-.4, .02)
    location <- c("bottomleft", "bottomright", "bottomright")
    x.axp.location <- list(c(-90, -52, -15, 0), c(-31, -15, 0, 15, 52, 59), c(-15, 0, 15, 52, 80))

    # pdf
    pdf.name <- paste(direct[1], ".pdf", sep = '')
    pdf(file = pdf.name, width = 12, height = 10)
    par(mfrow = c(3,3))

    for (m in 1:length(election)){
    elec <- get(file[(i - 1)*3 + m])
    loess.1 <- get(loess.slice.1[(i - 1)*3 + m])
    loess.2 <- get(loess.slice.2[(i - 1)*3 + m])

    elec$patet <- with(elec, widower.effect.est / control.turnout.est)
    loess.1$patet <- with(loess.1, widower.effect.est / control.turnout.est)
    loess.2$patet <- with(loess.2, widower.effect.est / control.turnout.est)

    mean.control.turnout <- round(mean(elec$control.turnout.est), 2)

    # case/control plot
    plot(elec$week, elec$case.turnout.est, main = election.name[m], xlab = x.title, ylab = "Turnout", xlim = x.axes[m], ylim = c(mean.control.turnout - .3, mean.control.turnout + .02), col = color[m], pch = 3, cex.lab = 1.15, bty = "n", xaxt = "n", yaxt = "n")
    points(elec$week, elec$control.turnout.est, bty = "n")
    legend(location[m], "(x,y)", c(filecontrols, filecases, match.ratio), pch = c(1, 3, 15), bty = "n", col = c("black", color[m], "black"), cex = 1.1)

    abline(v = 0, lty = 2)
    axis(side = 3, at = 0, tick = FALSE, padj = 1.5)
    axis(side = 1, at = x.axp.location[m])
    axis(side = 2, at = c(mean.control.turnout - .3, mean.control.turnout, 6))

    # loess plot
    color.loess <- rgb(col2rgb(color[m])[1,1], col2rgb(color[m])[2,1], col2rgb(color[m])[3,1], 175, maxColorValue=255)
    plot(elec$week, elec$patet, main = "", xlab = "", ylab = y.title, xlim = x.axes[m], ylim = y.axes, col = color.loess, pch = 16, cex.lab = 1.15, bty = "n", xaxt = "n")
    legend(location[m], c(filecases, match.ratio), pch = c(16, 15), bty = "n")

    abline(h = 0, lty = 2)
    abline(v = 0, lty = 2)
    axis(side = 1, at = x.axp.location[m])
    axis(side = 3, at = 0, tick = FALSE, padj = 1.5)

    season.dates <- c(40/7, 134/7, 223/7, 406/7, 499/7)
    seasons <- list(season.dates - 76, (season.dates - 21)[-2], season.dates)
    if (m == 2) {season.names <- c("<F", "<Sp","<F", "<Su")} else {season.names <- c("<F", "<Su", "<Sp","<F", "<Su")}
    year.dates <- c(-60/7, 305/7)
    years <- list(year.dates - 76, year.dates - 21, year.dates)
    year.names <- c("<2011", "<2010")

    axis(side = 3, at = seasons[m], tick = FALSE, padj = 1.5, labels = season.names)
    axis(side = 3, at = years[m], tick = FALSE, padj = 1.5, labels = year.names)

    if (m != 3) {lines(loess.smooth(loess.1$week, loess.1$patet, span = span.par[m][1]), lwd = 1.5, lty = 1)}
    if (m != 1) {lines(loess.smooth(loess.2$week, loess.2$patet, span = span.par[m][1]), lwd = 1.5, lty = 1)}

    # discontinuity regressions
    plot(elec$week, elec$patet, main = "", xlab = "", ylab = "", xlim = x.axes[m], ylim = y.axes, col = color.loess, bty = "n", xaxt = "n", pch = 16)
    legend(location[m], c(filecases, match.ratio), pch = c(3, 15), bty = "n")

    abline(h = 0, lty = 2)
    abline(v = 0, lty = 2)
    axis(side = 1, at = x.axp.location[m])
    axis(side = 3, at = 0, tick = FALSE, padj = 1.5)

    abline(v = 15, lty = 3)
    axis(side = 3, at = 15, tick = FALSE, padj = 1.5, labels = "15")
    abline(v = 52, lty = 3)
    axis(side = 3, at = 52, tick = FALSE, padj = 1.5, labels = "52")

    slope<-"slopes:"

    if (m != 1) {elec.15to52 <- elec[elec$week>=15&elec$week<=52,]
    elec.lm <- lm(patet ~ week, elec.15to52)
    summary.lm <- summary(elec.lm)
    lines(elec.15to52$week, elec.lm$fitted.values, lwd = 1.5)

    text(52, -.18, slope[1], cex = 1.15)
    slope[2]<-formatC(round(coef(elec.lm)[2], 5), digits = 1, format = "e")
    p.value <- paste("pvalue: ", round(coef(summary.lm)[2,4], 3), sep = "")
    text(33.5, -.235, slope[2], cex = 1.15)
    text(33.5, -.28, p.value, cex = 1)
    }

    if (m == 3) {
    elec.52to80 <- elec[elec$week>=52&elec$week<=80,]
    elec.lm <- lm(patet ~ week, elec.52to80)
    summary.lm <- summary(elec.lm)
    lines(elec.52to80$week, elec.lm$fitted.values, lwd = 1.5)

    slope[2]<-formatC(round(coef(elec.lm)[2], 5), digits = 1, format = "e")
    p.value <- paste("pvalue: ", round(coef(summary.lm)[2,4], 3), sep = "")
    text(67.5, -.215, slope[2], cex = 1.15)
    text(67.5, -.25, p.value, cex = 1)

    start <- mean(elec[elec$week>=15&elec$week<=25,]$widower.effect.est)
    end <- mean(elec[elec$week>=70&elec$week<=80,]$widower.effect.est)
    mtext(round(start, 2), side = 4, at = round(start, 2), las = 2, cex = .75, adj = .25)
    mtext(round(end, 2), side = 4, at = round(end, 2), las = 2, cex = .75, adj = .25)
    }
  }
    dev.off()
  }

## 2 Wide
## ~~~~~~~

## 2.1 Preamble
## =============
  rm(list = ls())

      option.name <- c("All Registered Voters")

      x.title <- "Weeks Since Spousal Death"
      y.title <- "Widower Effect"

      election <- c("ss9", "gp10", "gg10")

      election.name <- c("Special Statewide 2009", "Gubernatorial Primary 2010", "Gubernatorial General 2010")

## 2.2 Read files
## ===============
   file <- NULL
   agg.file <- NULL
   readfile <- NULL

   for (i in 1:length(option.name)) {
   for (j in 1:length(election)) {
      file[(i-1)*3 + j] <- paste(option.name[i], " ", election[j], sep = '') # filenames no ext.
      agg.file[(i - 1)*3 + j] <- paste(option.name[i], " - Aggregated ", election[j], sep = "")
  }}

     readfile <- paste("output/", file, ".csv", sep = '') # filenames + dir. & ext.
     file <- gsub(" ", ".", file)
     agg.readfile <- paste("output/", agg.file, ".csv", sep = "")
     agg.file <- gsub(" - ", ".", agg.file)
     agg.file <- gsub(" ", ".", agg.file)

     loess.slice <- paste(file, ".loess", sep = "")
     loess.slice.1 <- paste(file, ".loess.1", sep = "")
     loess.slice.2 <- paste(file, ".loess.2", sep = "")

      for (k in 1:length(readfile)){
       assign(file[k], read.csv(readfile[k], head = TRUE))
       assign(loess.slice.1[k], get(file[k])[get(file[k])$week<=0,])
       assign(loess.slice.2[k], get(file[k])[get(file[k])$week>=0,])

       assign(agg.file[k], read.csv(agg.readfile[k], head = TRUE))
     }

## 2.3 Plot
## =========
    for (i in 1:length(option.name)){
      direct <- paste('plots/Wide -- ', option.name[i], sep = '')

      # prepare legend variables
      filecases<-round(mean(get(file[(i - 1)*3 + 1])$cases.in.week.matched))
      filecases<-paste("# cases/week =", filecases[1])
      filecontrols<-round(mean(get(file[(i - 1)*3 + 1])$controltocase.ratio))
      filecontrols<-paste("controls:cases =", filecontrols[1])
      filecontrols<-paste(filecontrols[1], ":1", sep = '')
      match.ratio<-round(mean(get(file[(i - 1)*3 + 1])$match.ratio), 2)
      match.ratio <-paste("    match rate =", match.ratio[1])

      # prepare plot variables
      color <- c("goldenrod", "forestgreen", "purple4")
      color.highlight <- c(gray(.5), gray(.5), gray(.5))
      span.par <- list(c(.3, .3), c(.5, .4), c(.3, .3))
      x.axes <- list(c(-91, 4), c(-31, 59), c(-15, 80))
      y.axes <- c(-.25, .02)
      location <- c("bottomleft", "bottomright", "bottomright")
      x.axp.location <- c(-90, -52, -15, 0, 15, 52, 80)

      # pdf
      pdf.name <- paste(direct[1], ".pdf", sep = '')
      pdf(file = pdf.name, width = 8, height = 5, colormodel = "srgb")
      par(mfrow = c(1,1), mar = c(5, 4, 1, 2), family = "serif")

      for (m in 1:length(election)){
      elec <- get(file[(i - 1)*3 + m])
      agg.elec <- get(agg.file[(i - 1)*3 + m])
      loess.1 <- get(loess.slice.1[(i - 1)*3 + m])
      loess.2 <- get(loess.slice.2[(i - 1)*3 + m])

      elec$patet <- elec$widower.effect.est #/ elec$control.turnout.est
      loess.1$patet <- loess.1$widower.effect.est #/ loess.1$control.turnout.est
      loess.2$patet <- loess.2$widower.effect.est #/ loess.2$control.turnout.est
      agg.elec$patet <- agg.elec$widower.effect.est #/ agg.elec$control.turnout.est


      mean.control.turnout <- round(mean(elec$control.turnout.est), 2)

      # loess plot
      color.loess <- rgb(col2rgb(color[m])[1,1], col2rgb(color[m])[2,1], col2rgb(color[m])[3,1], 175, maxColorValue=255)
      if (m == 1) {plot(elec$week, elec$patet, main = "", xlab = x.title, ylab = "", xlim = c(-91, 80), ylim = y.axes, col = color.loess, pch = 16, cex = .8, cex.lab = .95, bty = "n", xaxt = "n", yaxt = "n", col.lab = gray(.2))


      avg.n90.to.n15 <- mean(agg.elec[agg.elec$week==1,]$patet)
      mtext(100*round(avg.n90.to.n15, 3), side = 2, at = avg.n90.to.n15, las = 2, cex = 1.15, adj = .6, col = color[m])

      segments(-91, avg.n90.to.n15, -52, avg.n90.to.n15, col = color[m], lwd = 2) # top-left long lines

                 }

      if (m != 1) {points(elec$week, elec$patet, col = color.loess, pch = 16, cex = .8)}

      points(elec$week[c(11,20)], elec$patet[c(11,20)], col = color.highlight[m], pch = 1, cex = 1.5, lwd = 2)
      points(elec$week[c(87,96)], elec$patet[c(87,96)], col = color.highlight[m], pch = 22, cex = 1.5, lwd = 2)
    }

      avg.15.to.80 <- mean(agg.elec[agg.elec$week==7,]$patet)
      mtext(100*round(avg.15.to.80, 3), side = 4, at = avg.15.to.80, las = 2, cex = 1.15, col = color[m], adj = .25)

      segments(52, avg.15.to.80, 80, avg.15.to.80, col = color[m], lwd = 2) # top-right long lines


  # add nadir average
      elec.ss9 <- get(file[(i - 1)*3 + 1])
      nadir.ss9 <- elec.ss9[elec.ss9$week>=0&elec.ss9$week<=2,]$widower.effect.est #/ elec.ss9[elec.ss9$week=>=0&elec.ss9$week<=3,]$control.turnout.est
      elec.gp10 <- get(file[(i - 1)*3 + 2])
      nadir.gp10 <- elec.gp10[elec.gp10$week>=0&elec.gp10$week<=2,]$widower.effect.est #/ elec.gp10[elec.gp10$week=>=0&elec.gp10$week<=3,]$control.turnout.est
      elec.gg10 <- get(file[(i - 1)*3 + 3])
      nadir.gg10 <- elec.gg10[elec.gg10$week>=0&elec.gg10$week<=2,]$widower.effect.est #/ elec.gg10[elec.gg10$week=>=0&elec.gg10$week<=3,]$control.turnout.est
      nadir <- (mean(nadir.ss9) + mean(nadir.gp10) + mean(nadir.gg10)) / 3

      text(10, nadir, 100*round(nadir, 3), cex = 1, col = gray(.3))

      abline(h = 0, lty = 2)
      abline(v = 0, lty = 2)
      axis(side = 2, at = 0, labels = "0.0", las = 2, tick = FALSE, hadj = .25)
      axis(side = 4, at = 0, labels = "0.0", las = 2, tick = FALSE, hadj = .75)

      mtext(y.title, side = 2, padj = -3.5, cex = 1, col = gray(.3), las = 0)
      axis(side = 1, at = x.axp.location, col.axis = gray(.2), tick = FALSE)
      legend("bottomleft", c("Special Statewide 2009", "Gubernatorial Primary 2010", "Gubernatorial General 2010"), pch = 16, col = color, bty = "n", cex = .8)

   dev.off()
    }


## 3 Gender comparison
## ~~~~~~~~~~~~~~~~~~~~

## 3.1 Preamble
## =============
  rm(list = ls())

      option.name <- c("Male Ego", "Female Ego")

      x.title <- "Weeks Since Spousal Death"
      y.title <- "Widower Effect"

      election <- c("ss9", "gp10", "gg10")

      election.name <- c("Special Statewide 2009", "Gubernatorial Primary 2010", "Gubernatorial General 2010")

## 3.2 Read files
## ===============
   file <- NULL
   agg.file <- NULL
   readfile <- NULL

   for (i in 1:length(option.name)) {
   for (j in 1:length(election)) {
      file[(i-1)*3 + j] <- paste(option.name[i], " ", election[j], sep = '') # filenames no ext.
      agg.file[(i - 1)*3 + j] <- paste(option.name[i], " - Aggregated ", election[j], sep = "")
  }}

     readfile <- paste("output/", file, ".csv", sep = '') # filenames + dir. & ext.
     file <- gsub(" ", ".", file)
     agg.readfile <- paste("output/", agg.file, ".csv", sep = "")
     agg.file <- gsub(" - ", ".", agg.file)
     agg.file <- gsub(" ", ".", agg.file)

     loess.slice <- paste(file, ".loess", sep = "")
     loess.slice.1 <- paste(file, ".loess.1", sep = "")
     loess.slice.2 <- paste(file, ".loess.2", sep = "")

      for (k in 1:length(readfile)){
       assign(file[k], read.csv(readfile[k], head = TRUE))
       assign(loess.slice.1[k], get(file[k])[get(file[k])$week<=0,])
       assign(loess.slice.2[k], get(file[k])[get(file[k])$week>=0,])

       assign(agg.file[k], read.csv(agg.readfile[k], head = TRUE))
     }

## 3.3 Plot
## =========
    pdf(width = 7.5, height = 5, file = "plots/Comparison -- Gender.pdf", colormodel = "srgb")
    par(family = "serif", fg = gray(.1), mar = c(4.5, 4, 1.25, 3)) # set font and colors

     # dummy points (white)
     with(get(file[k]), plot(week, widower.effect.est / control.turnout.est, xlab = x.title, ylab = "", xlim = c(-91, 80), ylim = c(-.25, .05), col = "white", xaxp = c(-90, 80, 17), yaxp = c(-.45, .05, 10) , pch = 3, cex.lab = 1, bty = "n", xaxt = "n", yaxt = "n", col.lab = gray(.2), cex.main = 1.3))

     # add lines
     abline(h = 0, lty = 2, col = gray(.1))  # 0 for no change in turnout
     abline(v = 0, lty = 2, col = gray(.1)) # 0 for week 0
     axis(side = 1, at = c(-90, -52, -15, 0, 15, 52, 80), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))    # x axis labels
     mtext("0.0", side = 4, at = 0,  las = 2)

     # add subtitle
     mtext(y.title, side = 2, padj = -3.75, cex = 1.15, col = gray(.3), las = 0)

     # label above means
     mtext(" weeks
   -91 to -52", side = 3, at = -98, padj = 1.1, cex = .9, col = gray(.3))
     mtext("  weeks
    52 to 80", side = 3, at = 90, padj = 1.1, cex = .9, col = gray(.3))

     # mention that these lines are from two elections
     mtext("LHS -- SS 2009", side = 1, at = -70, padj = 4, col = "goldenrod")
     mtext("RHS -- GG 2010", side = 1, at = 60, padj = 4, col = "purple3")

    # line-specifc variables
    color <- c("steelblue4", "plum3")
    line.type <- c(1, 1)    # more, same, fewer
    line.width <- c(2.25, 2.25)
    option.name <- gsub("Ego", "Subject", option.name)

    # draw the legend
    legend("bottomleft", "(x,y)", option.name, lty = line.type, cex = 1.1, text.col = color, bty = "n", lwd = line.width, col = color)

    # plot loess lines
    for (m in 1:length(option.name)){

      elec <- get(loess.slice.1[(m-1)*3 + 1])
      agg.elec <- get(agg.file[(m-1)*3 + 1])
      elec$patet <- elec$widower.effect.est #/ elec$control.turnout.est
      agg.elec$patet <- agg.elec$widower.effect.est #/ agg.elec$control.turnout.est
      lines(loess.smooth(elec$week, elec$patet, span = .3), col = color[m], lty = line.type[m], lwd = line.width[m])
      avg.n90.to.n52 <- mean(agg.elec[agg.elec$week==1,]$patet)

      # add weeks -91 to -52 means
      mtext(sprintf("%.1f", 100*round(avg.n90.to.n52, 2)), side = 2, at = round(avg.n90.to.n52, 2), las = 2, cex = .9, adj = .6, col = color[m])

      elec <- get(loess.slice.2[(m-1)*3 + 3])
      agg.elec <- get(agg.file[(m-1)*3 + 3])
      elec$patet <- elec$widower.effect.est #/ elec$control.turnout.est
      agg.elec$patet <- agg.elec$widower.effect.est #/ agg.elec$control.turnout.est
      lines(loess.smooth(elec$week, elec$patet, span = .3), col = color[m], lty = line.type[m], lwd = line.width[m])
      avg.52.to.80 <- mean(agg.elec[agg.elec$week==7,]$patet)

      # add weeks 52 to 80 means
      mtext(sprintf("%.1f", 100*round(avg.52.to.80, 3)), side = 4, at = round(avg.52.to.80, 3), las = 2, cex = .9, col = color[m], adj = .25)

      #  # these segments indicate the averaged weeks
    segments(-91, .02, -52, .02, col = gray(.4)) # top-left long lines
    segments(52, .02, 80, .02, col = gray(.4)) # top-right long lines

    }

  dev.off()

## 4 Ages
## ~~~~~~~

## 4.1 Preamble
## =============
  rm(list = ls())

      option.name <- c("Age - Same Voting History", "Age - More Votes than Alter", "Age - Fewer Votes than Alter")

      x.title <- "Ego Age"
      y.title <- "Widower Effect"
      y.sub.title <- "aggregated for weeks 15 through 80"
      election <- c("gg10", "gp10", "ss9")

## 4.2 Read files
## ===============
   file <- NULL
   readfile <- NULL
   patet <- NULL

   for (i in 1:length(option.name)) {
      ## file[i] <- paste(option.name[1], " ", election[i], sep = '') # filenames no ext.
      file[i] <- paste(option.name[i], " ", election[1], sep = '') # filenames no ext.
  }

     readfile <- paste("output/", file, ".csv", sep = '') # filenames + dir. & ext.
     file <- gsub(" ", ".", file)

      for (k in 1:length(readfile)){
       assign(file[k], read.csv(readfile[k], head = TRUE))
     }

## 4.3 Same voting history plot
## =============================
      pdf(width = 7, height = 5, file = "plots/By -- age.pdf", colormodel = "srgb")
      par(family = "serif", fg = gray(.1), mar = c(4.5, 6, 1.25, 3)) # set font and colors

       # dummy points (white)
       with(get(file[k]), plot(ego.age.avg, 100*widower.effect.est / control.turnout.est, xlab = x.title, ylab = "", xlim = c(45, 95), ylim = c(-25, 3), col = "white", xaxp = c(-90, 80, 17), yaxp = c(-45, 5, 10) , pch = 3, cex.lab = 1, bty = "n", xaxt = "n", col.lab = gray(.2), cex.main = 1.3))

       # add lines
       abline(h = 0, lty = 2, col = gray(.1))  # 0 for no change in turnout
       axis(side = 1, at = c(45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))    # x axis labels
       mtext("0.0", side = 4, at = 0,  las = 2)

       # add y title and subtitle
       mtext(y.sub.title, side = 2, padj = -5, cex = 1, col = gray(.2), las = 0)
       mtext(y.title, side = 2, padj = -6, cex = 1.15, col = gray(.3), las = 0)

      # line-specifc variables
      color <- c("purple4", "mediumpurple", "purple4")
      line.type <- c(1, 2, 6)    # same, more, fewer
      line.width <- c(2.25, 2.25, 2.25)

      # draw the legend
      option.name.plot <- gsub("Votes than Alter", "votes than deceased", option.name)
      option.name.plot <- gsub("Same Voting History", "Same voting history", option.name.plot)
      option.name.plot <- gsub("Age - ", "", option.name.plot)
      legend("bottomright", "(x,y)", option.name.plot[m], lty = line.type, cex = 1.1, text.col = color, bty = "n", lwd = line.width, col = color[m])

      # plot loess lines
      for (m in 1:1){

        elec <- get(file[m])
        elec <- elec[elec$ego.age.avg>=45,]
        elec$patet <- 100 * elec$widower.effect.est #/ elec$control.turnout.est
        points(elec$ego.age.avg, 100 * elec$widower.effect.est, col = color[m], pch = 3)
  #      points(elec$ego.age.avg, elec$control.turnout.est, col = "black", pch = 3)
        lines(loess.smooth(elec$ego.age.avg, elec$patet, span = .4), col = color[m], lty = line.type[m], lwd = line.width[m])

      }

    dev.off()

## 4.4 Voting history comparison plot
## ===================================
    pdf(width = 7, height = 5, file = "plots/Comparison -- Age.pdf", colormodel = "srgb")
    par(family = "serif", fg = gray(.1), mar = c(4.5, 6, 1.25, 3)) # set font and colors

     # dummy points (white)
     with(get(file[k]), plot(ego.age.avg, 100 * widower.effect.est / control.turnout.est, xlab = x.title, ylab = "", xlim = c(45, 95), ylim = c(-25, 0), col = "white", xaxp = c(-90, 80, 17), yaxp = c(-45, 5, 10) , pch = 3, cex.lab = 1, bty = "n", xaxt = "n", col.lab = gray(.2), cex.main = 1.3))

     # add lines
     abline(h = 0, lty = 2, col = gray(.1))  # 0 for no change in turnout
     axis(side = 1, at = c(45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))    # x axis labels
     mtext("0.0", side = 4, at = 0,  las = 2)

     # add y title and subtitle
     mtext(y.sub.title, side = 2, padj = -5, cex = 1, col = gray(.2), las = 0)
     mtext(y.title, side = 2, padj = -6, cex = 1.15, col = gray(.3), las = 0)

    # line-specifc variables
    color <- c("purple4", "mediumpurple", "purple4")
    line.type <- c(6, 1, 2)    # more, same, fewer
    line.width <- c(2.25, 2.25, 2.25)

    # draw the legend
    option.name.plot <- gsub("Votes than Alter", "votes than deceased", option.name)
    option.name.plot <- gsub("Same Voting History", "Same voting history", option.name.plot)
    option.name.plot <- gsub("Age - ", "", option.name.plot)
    legend("bottomright", "(x,y)", option.name.plot[2:3], lty = line.type[2:3], cex = 1.1, text.col = color[2:3], bty = "n", lwd = line.width[2:3], col = color[2:3])

    # plot loess lines
    for (m in 2:length(option.name)){

      elec <- get(file[m])
      elec <- elec[elec$ego.age.avg>=45,]
      elec$patet <- 100 * elec$widower.effect.est #/ elec$control.turnout.est
      points(elec$ego.age.avg, elec$patet, col = color[m], pch = 3)
      lines(loess.smooth(elec$ego.age.avg, elec$patet, span = .4), col = color[m], lty = line.type[m], lwd = line.width[m])

    }

  dev.off()

## 4.5 Age and turnout
## ====================

## 4.5.1 Treatment and age
## ------------------------
  agg.sample <- read.csv("output/aggregated_overall_sample_turnout.csv")
  agg.sample.same <- read.csv("output/aggregated_overall_same_sample_turnout.csv")

  agg.sample <- agg.sample[agg.sample$age.gg10.e>=24&agg.sample$age.gg10.e<=99,]
  agg.sample.same <- agg.sample.same[agg.sample.same$age.gg10.e>=24&agg.sample.same$age.gg10.e<=99,]

  pdf(width = 6, height = 8, file = "plots/Age and past turnout.pdf", colormodel = "srgb")
  par(mfrow = c(3, 2))

  attach(agg.sample.same, warn.conflicts = FALSE)
  plot(age.gg10.e, gg10.e, pch = 16, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", main = "Gubernatorial General 2010", cex.main = .9, cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(gg10.e), 2), round(max(gg10.e), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", "Full spouse sample", pch = 16, bty = "n")

  attach(get(file[1]), warn.conflicts = FALSE)
  color <- "purple4"
  color <- rgb(col2rgb(color)[1,1], col2rgb(color)[2,1], col2rgb(color)[3,1], 175, maxColorValue=255)
  plot(ego.age.avg, case.turnout.est, pch = 16, col = color, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", cex = .8)
  points(ego.age.avg, control.turnout.est, pch = 16, col = "black", xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty= "n", cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(control.turnout.est), 2), round(max(control.turnout.est), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", c("Matched spouses", "Widow(er)s"), pch = 16, col = c("black", color), bty = "n")



  attach(agg.sample.same, warn.conflicts = FALSE)
  plot(age.gg10.e, gp10.e, pch = 16, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", main = "Gubernatorial Primary 2010", cex.main = .9, cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(gp10.e), 2), round(max(gp10.e), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", "Full spouse sample", pch = 16, bty = "n")

  attach(get(file[2]), warn.conflicts = FALSE)
  color <- "forestgreen"
  color <- rgb(col2rgb(color)[1,1], col2rgb(color)[2,1], col2rgb(color)[3,1], 175, maxColorValue=255)
  plot(age.avg, case.turnout.est, pch = 16, col = color, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", cex = .8)
  points(age.avg, control.turnout.est, pch = 16, col = "black", xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(control.turnout.est), 2), round(max(control.turnout.est), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", c("Matched spouses", "Widow(er)s"), pch = 16, col = c("black", color), bty = "n")



  attach(agg.sample.same, warn.conflicts = FALSE)
  plot(age.gg10.e, ss9.e, pch = 16, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", main = "Special Statewide 2009", cex.main = .9, cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(ss9.e), 2), round(max(ss9.e), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", "Full spouse sample", pch = 16, bty = "n")

  attach(get(file[3]), warn.conflicts = FALSE)
  color <- "goldenrod"
  color <- rgb(col2rgb(color)[1,1], col2rgb(color)[2,1], col2rgb(color)[3,1], 175, maxColorValue=255)
  plot(age.avg, case.turnout.est, pch = 16, col = color, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", cex = .8)
  points(age.avg, control.turnout.est, pch = 16, col = "black", xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(control.turnout.est), 2), round(max(control.turnout.est), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", c("Matched spouses", "Caregivers"), pch = 16, col = c("black", color), bty = "n")

  dev.off()

## 4.5.2 Control and age
## ----------------------
  agg.sample <- read.csv("output/aggregated_overall_sample_turnout.csv")
  agg.sample.same <- read.csv("output/aggregated_overall_same_sample_turnout.csv")

  agg.sample <- agg.sample[agg.sample$age.gg10.e>=24&agg.sample$age.gg10.e<=99,]
  agg.sample.same <- agg.sample.same[agg.sample.same$age.gg10.e>=24&agg.sample.same$age.gg10.e<=99,]

  pdf(width = 6, height = 8 * 4/3, file = "plots/Age and past turnout - control elections.pdf", colormodel = "srgb")
  par(mfrow = c(4, 2))

  attach(agg.sample.same, warn.conflicts = FALSE)
  plot(age.gg10.e - 4, gg6.e, pch = 16, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", main = "Gubernatorial General 2006", cex.main = .9, cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(gg6.e), 2), round(max(gg6.e), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", "Full spouse sample", pch = 16, bty = "n")

  attach(get(file[1]), warn.conflicts = FALSE)
  color <- "purple4"
  color <- rgb(col2rgb(color)[1,1], col2rgb(color)[2,1], col2rgb(color)[3,1], 175, maxColorValue=255)
  plot(ego.age.avg - 4, turnout.gg6, pch = 16, col = color, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(turnout.gg6), 2), round(max(turnout.gg6), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", c("Exactly matched turnout"), pch = 16, col = c(color), bty = "n")



  attach(agg.sample.same, warn.conflicts = FALSE)
  plot(age.gg10.e - 4, gp6.e, pch = 16, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", main = "Gubernatorial Primary 2006", cex.main = .9, cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(gp6.e), 2), round(max(gp6.e), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", "Full spouse sample", pch = 16, bty = "n")

  attach(get(file[2]), warn.conflicts = FALSE)
  color <- "forestgreen"
  color <- rgb(col2rgb(color)[1,1], col2rgb(color)[2,1], col2rgb(color)[3,1], 175, maxColorValue=255)
  plot(age.avg - 4, turnout.gp6, pch = 16, col = color, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(turnout.gp6), 2), round(max(turnout.gp6), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", c("Exactly matched turnout"), pch = 16, col = c(color), bty = "n")


  attach(agg.sample.same, warn.conflicts = FALSE)
  plot(age.gg10.e -  4, ss5.e, pch = 16, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", main = "Special Statewide 2005", cex.main = .9, cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(ss5.e), 2), round(max(ss5.e), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", "Full spouse sample", pch = 16, bty = "n")

  attach(get(file[3]), warn.conflicts = FALSE)
  color <- "goldenrod"
  color <- rgb(col2rgb(color)[1,1], col2rgb(color)[2,1], col2rgb(color)[3,1], 175, maxColorValue=255)
  plot(age.avg - 4, turnout.ss5, pch = 16, col = color, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(turnout.ss5), 2), round(max(turnout.ss5), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", c("Exactly matched turnout"), pch = 16, col = c(color), bty = "n")

  attach(agg.sample.same, warn.conflicts = FALSE)
  plot(age.gg10.e -  4, pg4.e, pch = 16, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", main = "Presidential General 2004", cex.main = .9, cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(pg4.e), 2), round(max(pg4.e), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", "Full spouse sample", pch = 16, bty = "n")

  attach(get(file[3]), warn.conflicts = FALSE)
  color <- "deeppink"
  color <- rgb(col2rgb(color)[1,1], col2rgb(color)[2,1], col2rgb(color)[3,1], 175, maxColorValue=255)
  plot(age.avg - 4, turnout.pg4, pch = 16, col = color, xlab = "Age", ylab = "Turnout", ylim = c(0, 1), xlim = c(24, 99), bty = "n", xaxt = "n", yaxt = "n", cex = .8)
  axis(side = 1, at = c(25, 45, 55, 65, 75, 85, 95), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))
  axis(side = 2, at = c(round(min(turnout.pg4), 2), round(max(turnout.pg4), 2)), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), lwd = 0, lwd.ticks = 1)
  axis(side = 2, at = c(0, 1), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1), labels = FALSE)
  legend("bottomright", c("Exactly matched turnout"), pch = 16, col = c(color), bty = "n")

  dev.off()

## 5 Voting History Comparison
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~

## 5.1 Preamble
## =============
  rm(list = ls())

  plot.name <- "Voting History Comparison"

  option.name <- c("More Votes than Alter", "Same Voting History", "Fewer Votes than Alter")

  x.title <- "Weeks Since Spousal Death"
  y.title <- "Widower Effect"

  election <- c("ss9", "gp10", "gg10")

## 5.2 Read files
## ===============
   file <- NULL
   agg.file <- NULL
   readfile <- NULL

   for (i in 1:length(option.name)) {
   for (j in 1:length(election)) {
      file[(i-1)*3 + j] <- paste(option.name[i], " ", election[j], sep = '') # filenames no ext.
      agg.file[(i - 1)*3 + j] <- paste(option.name[i], " - Aggregated ", election[j], sep = "")
  }}

     readfile <- paste("output/", file, ".csv", sep = '') # filenames + dir. & ext.
     file <- gsub(" ", ".", file)
     agg.readfile <- paste("output/", agg.file, ".csv", sep = "")
     agg.file <- gsub(" - ", ".", agg.file)
     agg.file <- gsub(" ", ".", agg.file)

     loess.slice <- paste(file, ".loess", sep = "")
     loess.slice.1 <- paste(file, ".loess.1", sep = "")
     loess.slice.2 <- paste(file, ".loess.2", sep = "")

      for (k in 1:length(readfile)){
       assign(file[k], read.csv(readfile[k], head = TRUE))
       assign(loess.slice.1[k], get(file[k])[get(file[k])$week<=0,])
       assign(loess.slice.2[k], get(file[k])[get(file[k])$week>=0,])

       assign(agg.file[k], read.csv(agg.readfile[k], head = TRUE))
     }

## 5.3 Plot
## =========
    pdf(width = 7.5, height = 5, file = "plots/Comparison -- Voting.pdf", colormodel = "srgb")
    par(family = "serif", fg = gray(.1), mar = c(4.5, 4, 1.25, 3)) # set font and colors

     # dummy points (white)
     elec <- get(file[1])
     plot(elec$week, elec$widower.effect.est, xlab = x.title, ylab = "", xlim = c(-91, 80), ylim = c(-.3, .05), col = "white", xaxp = c(-90, 80, 17), yaxp = c(-.45, .05, 10) , pch = 3, cex.lab = 1, bty = "n", xaxt = "n", yaxt = "n", col.lab = gray(.2), cex.main = 1.3)

     # add lines
     abline(h = 0, lty = 2, col = gray(.1))  # 0 for no change in turnout
     abline(v = 0, lty = 2, col = gray(.1)) # 0 for week 0
     axis(side = 1, at = c(-90, -52, -15, 0, 15, 52, 80), col = gray(.1), col.axis = gray(.2), col.ticks = gray(.1))    # x axis labels
     mtext("0.0", side = 4, at = 0,  las = 2)

     # add subtitle
     mtext(y.title, side = 2, padj = -3.75, cex = 1.15, col = gray(.3), las = 0)

     # label above means
     mtext(" weeks
   -91 to -52", side = 3, at = -98, padj = 1.1, cex = .9, col = gray(.3))
     mtext("  weeks
    52 to 80", side = 3, at = 90, padj = 1.1, cex = .9, col = gray(.3))

     # mention that these lines are from two elections
     mtext("LHS -- SS 2009", side = 1, at = -70, padj = 4, col = "goldenrod")
     mtext("RHS -- GG 2010", side = 1, at = 60, padj = 4, col = "purple4")

    # line-specifc variables
    color <- c(gray(.6), gray(.4), gray(.2))    # it was colorful once
    color.1 <- c("mediumpurple", gray(.45), "purple4")
    color.2 <- c("goldenrod2", gray(.6), "goldenrod4")
    line.type <- c(2, 1, 6)    # more, same, fewer
    line.width <- c(2.25, 2.25, 2.25)

    # draw the legend
    option.name.plot <- gsub("More Votes than Alter", "Deceased spouse voted less", option.name)
    option.name.plot <- gsub("Fewer Votes than Alter", "Deceased spouse voted more", option.name.plot)
    option.name.plot <- gsub("Same Voting History", "Same voting history", option.name.plot)
    legend("bottomleft", "(x,y)", option.name.plot, lty = line.type, cex = 1.1, text.col = color, bty = "n", lwd = line.width, col = color)

    # plot loess lines
    for (m in 1:length(option.name)){

      elec <- get(loess.slice.1[(m-1)*3 + 1])
      min.elec.week <- min(elec$week)
      max.elec.week <- max(elec$week)
      original.elec.patet <- elec$widower.effect.est
      elec$week <- cut(elec$week, breaks = seq(min.elec.week, max.elec.week, 4), right = F)
      elec <- aggregate(widower.effect.est ~ week, data = elec, FUN = mean)
      agg.elec <- get(agg.file[(m-1)*3 + 1])
      elec$patet <- elec$widower.effect.est #/ elec$control.turnout.est
      agg.elec$patet <- agg.elec$widower.effect.est #/ agg.elec$control.turnout.est
      points(seq(min.elec.week + 2, max.elec.week - 2, 4), elec$patet, col = color.2[m], lty = line.type[m], lwd = line.width[m], pch = 3)
      lines(loess.smooth(seq(min.elec.week, max.elec.week, 1), original.elec.patet, span = .3), col = color.2[m], lty = line.type[m], lwd = line.width[m])
      avg.n90.to.n52 <- mean(agg.elec[agg.elec$week==1,]$patet)


      # add weeks -91 to -52 means
      mtext(sprintf("%.1f", 100*round(avg.n90.to.n52, 3)), side = 2, at = round(avg.n90.to.n52, 2), las = 2, cex = .9, adj = .6, col = color.2[m])

      elec <- get(loess.slice.2[(m-1)*3 + 3])
      min.elec.week <- min(elec$week)
      max.elec.week <- max(elec$week)
      original.elec.patet <- elec$widower.effect.est
      elec$week <- cut(elec$week, breaks = seq(min.elec.week, max.elec.week, 4), right = F)
      elec <- aggregate(widower.effect.est ~ week, data = elec, FUN = mean)
      agg.elec <- get(agg.file[(m-1)*3 + 3])
      elec$patet <- elec$widower.effect.est #/ elec$control.turnout.est
      agg.elec$patet <- agg.elec$widower.effect.est #/ agg.elec$control.turnout.est
      points(seq(min.elec.week + 2, max.elec.week - 2, 4), elec$patet, col = color.1[m], lty = line.type[m], lwd = line.width[m], pch = 3)
      lines(loess.smooth(seq(min.elec.week, max.elec.week, 1), original.elec.patet, span = .3), col = color.1[m], lty = line.type[m], lwd = line.width[m])
      avg.52.to.80 <- mean(agg.elec[agg.elec$week==7,]$patet)

      # add weeks 52 to 80 means
      mtext(sprintf("%.1f", 100*round(avg.52.to.80, 3)), side = 4, at = round(avg.52.to.80, 2), las = 2, cex = .9, col = color.1[m], adj = .25)

      #  # these segments indicate the averaged weeks
    segments(-91, .02, -52, .02, col = gray(.4)) # top-left long lines
    segments(52, .02, 80, .02, col = gray(.4)) # top-right long lines

    }

  dev.off()



## 6 Recovery Plot
## ~~~~~~~~~~~~~~~~

## 6.1 Preamble
## =============
  rm(list = ls())

  #  option.name <- list(c("Alone", "Not Alone"), c("Rural Zip Code", "Urban Zip Code"), c("Ego Under 65", "Ego 65 to 79", "Ego 80 or Over"), c("Per Capita Income Under 35000", "Per Capita Income 35000 or Over"), c("More Votes than Alter", "Fewer Votes than Alter"), c("Not Alone", "Movers Only"))
  #      option.title <- c("Alone vs Not Alone", "Population Density", "Age", "Per Capita Income", "More vs Fewer", "Not Alone vs Movers Only")
    option.name <- list(c("Alone", "Not Alone"), c("Ego Under 65", "Ego 65 to 79", "Ego 80 or Over"))
        option.title <- c("Alone vs Not Alone", "Age")

        x.title <- "Weeks Since Spousal Death"
        y.title <- "Proportional Widower Effect"

        election <- c("ss9", "gp10", "gg10")

        election.name <- c("Special Statewide 2009", "Gubernatorial Primary 2010", "Gubernatorial General 2010")


## 6.2 Read files
## ===============
  file <- NULL
         readfile <- NULL
        move.file <- NULL
          for (i in 1:length(unlist(option.name))) {

         for (j in 1:length(election)) {
            file[(i - 1)*3 + j] <- paste(unlist(option.name)[i], " ", election[j], sep = '') # filenames no ext.
        }
            move.file[i] <- paste(unlist(option.name)[i], " - Movers", sep = '')
        }
           readfile <- paste("output/", file, ".csv", sep = '') # filenames + dir. & ext.
           move.readfile <- paste("output/", move.file, " gg10.csv", sep = '')
           file <- gsub(" ", ".", file)
           move.file <- gsub(" - ", ".", move.file); move.file <- gsub(" ", ".", move.file)

           loess.slice.1 <- paste(file, ".loess.1", sep = "")
           loess.slice.2 <- paste(file, ".loess.2", sep = "")

            for (k in 1:length(readfile)){
             assign(file[k], read.csv(readfile[k], head = TRUE))
             assign(loess.slice.1[k], get(file[k])[get(file[k])$week<=0,])
             assign(loess.slice.2[k], get(file[k])[get(file[k])$week>=0,])
           }
            for (k in 1:length(move.readfile)){
             assign(move.file[k], read.csv(move.readfile[k], head = TRUE))
             }

## 6.3 Plot
## =========
      direct <- paste('plots', "Recovery -- ", sep = '/')

      # prepare plot variables
      color <- "purple4"
      color <- rgb(col2rgb(color)[1,1], col2rgb(color)[2,1], col2rgb(color)[3,1], 50, maxColorValue=255)
      color.anniversary <- rgb(col2rgb(color)[1,1], col2rgb(color)[2,1], col2rgb(color)[3,1], 150, maxColorValue=255)
      x.axes <- c(15, 80)
      y.axes <- c(-.25, 0)
      text.location <- list()
      x.axp.location <- c(15, 52, 80)

    for (j in 1:length(option.title)) {

      # pdf
      pdf.name <- paste(direct[1], option.title[j], ".pdf", sep = '')
  if (length(option.name[j])==2){
      pdf(file = pdf.name, width = 7, height = 3, family = "serif")
      par(mfrow = c(1,2), mar = c(0, 2, 1, 2))
  } else {
      pdf(file = pdf.name, width = 7.5, height = 2.25, family = "serif")
      par(mfrow = c(1,3), mar = c(0, 2, 1, 2))}

    for (i in 1:length(option.name[j])){

      for (m in 3:length(election)){

      file.name <- paste("^", gsub(" ", ".", option.name[j][i]), ".gg10", sep = "")
      file.number <- grep(file.name, file)
      elec <- get(file[file.number])
      elec <- elec[elec$week>=15,]
      elec$patet <- elec$widower.effect.est #/ elec$control.turnout.est
          y.title <- "Widower Effect"

      # discontinuity regressions
      if (i <=2 ) {plot(elec$week, elec$patet, main = "", xlab = "", ylab = "", xlim = x.axes, ylim = y.axes, col = color, bty = "n", xaxt = "n", yaxt = "n", pch = 16, cex = .75)}
      if (i > 2 ) {plot(elec$week, elec$patet, main = "", xlab = "", ylab = "", xlim = x.axes, ylim = y.axes, col = color, bty = "n", xaxt = "n", yaxt = "n", pch = 16, cex = .75)}
      anniversary <- elec[elec$week==52,]
      points(anniversary$week, anniversary$patet, col = color.anniversary, pch = 16, cex = 1)
    #  legend(location, c(filecases, match.ratio), pch = c(3, 15), bty = "n")

      text(52, -.01, option.name[j][i], cex = .9, col = gray(.1))
   #   if (i > 2 ) {text(52, -.045, option.name[j][i], cex = .9, col = gray(.1))}

      segments(52, (y.axes[2] - .025), 52, (y.axes[1] + .03), lty = 3)
      axis(side = 1, at = 52, col = gray(.1), labels = "One year anniversary", padj = -5, cex.axis = .75, tick = FALSE)
      axis(side = 1, at = c(15, 52, 80), col = gray(.2), labels = c("15", "52", "80"), padj = -9, cex.axis = .6, tick = FALSE)

      slope<-"slopes:"

      elec.15to52 <- elec[elec$week>=15&elec$week<=52,]
      elec.lm <- lm(patet ~ week, elec.15to52)
      summary.lm <- summary(elec.lm)

      out.of.state.ratio.1 <- .16 / (1 - .16)
      out.of.state.ratio.2 <- .20 / (1 - .20)
      out.of.state.1 <- out.of.state.ratio.1*mean(elec.15to52$control.turnout.est)
      out.of.state.2 <- out.of.state.ratio.2*mean(elec.15to52$control.turnout.est)

      move.elec <- gsub(".gg10", ".Movers", file.name)
      move.elec.number <- grep(move.elec, move.file)
      move.elec.15 <- loess.smooth(get(move.file[move.elec.number ])$week, get(move.file[move.elec.number ])$widower.effect.est)$y[16]
      move.elec.52 <- loess.smooth(get(move.file[move.elec.number ])$week, get(move.file[move.elec.number ])$widower.effect.est)$y[36]
      move.elec.80 <- loess.smooth(get(move.file[move.elec.number ])$week, get(move.file[move.elec.number ])$widower.effect.est)$y[50]

     mover.adj <- seq(move.elec.15*out.of.state.1, move.elec.52*out.of.state.1, (move.elec.52*out.of.state.1 - move.elec.15*out.of.state.1) / (52 - 15))
     mover.adj.lm <- mover.adj + elec.lm$fitted.values
     lines(elec.15to52$week, mover.adj.lm, lwd = 1.5, lty = 2)
     lines(elec.15to52$week, elec.lm$fitted.values, lwd = 1.5)

    #  text(52, -.15, slope[1], cex = 1.15)
      slope[2]<- 100 * round(coef(elec.lm)[2], 5)
      p.value <- paste("pvalue: ", round(coef(summary.lm)[2,4], 3), sep = "")
      text(33.5, (y.axes[1] + .0675), slope[2], cex = .7, col = gray(.1))
      text(33.5, (y.axes[1] + .05), p.value, cex = .65, col = gray(.1))

      elec.52to80 <- elec[elec$week>=52&elec$week<=80,]
      elec.lm <- lm(patet ~ week, elec.52to80)
      summary.lm <- summary(elec.lm)

      mover.adj <- seq(move.elec.52*out.of.state.1, move.elec.80*out.of.state.2, (move.elec.80*out.of.state.2 - move.elec.52*out.of.state.1) / (80 - 52))
      mover.adj.lm <- mover.adj + elec.lm$fitted.values

      lines(elec.52to80$week, mover.adj.lm, lwd = 1.5, lty = 2)
      lines(elec.52to80$week, elec.lm$fitted.values, lwd = 1.5)

      slope[2]<-100 * round(coef(elec.lm)[2], 5)
      p.value <- paste("pvalue: ", round(coef(summary.lm)[2,4], 3), sep = "")
      text(67.5, (y.axes[1] + .0675), slope[2], cex = .7, col = gray(.1))
      text(67.5, (y.axes[1] + .05), p.value, cex = .65, col = gray(.1))

      start <- mean(elec[elec$week>=15&elec$week<=25,]$widower.effect.est)
      end <- mean(elec[elec$week>=70&elec$week<=80,]$widower.effect.est)
      mtext(sprintf("%.1f", 100*round(start, 3)), side = 2, at = round(start, 2), las = 2, cex = .75, adj = .9)
      mtext(sprintf("%.1f", 100*round(end, 3)), side = 4, at = round(end, 2), las = 2, cex = .75, adj = .25)
      end.adj <- mover.adj.lm[length(mover.adj.lm)] - elec.lm$fitted.values[length(elec.lm$fitted.values)]
      mtext(sprintf("%.1f", 100*round(end + end.adj, 3)), side = 4, at = round(end + end.adj, 2), las = 2, cex = .75, adj = .25)
    }
    }
    dev.off()
  }

## 7 Aggregate numbers
## ~~~~~~~~~~~~~~~~~~~~

## 7.1 Preamble
## =============
  rm(list = ls())

    option.name <- c("All Registered Voters", "Female Ego", "Male Ego", "Ego Under 65", "Ego 65 to 79", "Ego 80 or Over", "Female Ego Under 65", "Male Ego Under 65", "Female Ego 65 to 79", "Male Ego 65 to 79", "Female Ego 80 or Over", "Male Ego 80 or Over", "Female Ego - Opposite Sex Only", "Male Ego - Opposite Sex Only", "More Votes than Spouse", "Fewer Votes than Spouse", "Same Voting History", "Same Voting History - No Abstains", "Same Voting History - One Abstain", "Same Voting History - Two Abstains", "Same Voting History - Three Abstains", "Alone", "Not Alone", "Rural Zip Code", "Urban Zip Code", "Per Capita Income Under 35000", "Per Capita Income 35000 or Over", "California born", "Not California born")

    column.names <- c("weeks -91 to 52", "weeks -52 to -15", "weeks -15 to 0", "weeks 0 to 15", "weeks 15 to 52", "weeks 52 to 80")

    election <- c("ss9", "gp10", "gg10")

    election.name <- c("Special Statewide 2009", "Gubernatorial Primary 2010", "Gubernatorial General 2010")

## 7.2 Read files
## ===============
   file <- NULL
   readfile <- NULL

    for (i in 1:length(option.name)) {

   for (j in 1:length(election)) {
      file[(i - 1)*3 + j] <- paste(option.name[i], " - Aggregated ", election[j], sep = '') # filenames no ext.
  }
  }
     readfile <- paste("output/", file, ".csv", sep = '') # filenames + dir. & ext.
     file <- gsub(" ", ".", file)
     file <- gsub("-.", "", file)

      for (k in 1:length(readfile)){
       assign(file[k], read.csv(readfile[k], head = TRUE))
     }

## 7.3 Aggregate
## ==============
  elec.data <- NULL
  table <- NULL
  neg.91.to.neg.52 <- NULL
  neg.52.to.neg.15 <- NULL
  neg.15.to.0 <- NULL
  x.0.to.15 <- NULL
  x.15.to.52 <- NULL
  x.52.to.80 <- NULL
  blank <- rbind("", "", "")
  i <- NULL
  m <- NULL

  for (i in 1:length(option.name)){

      n <- (i - 1)*3

      m <- 1
      elec.data <- get(file[n + m])

      neg.91.to.neg.52 <- paste("n = ", elec.data[elec.data$week==1,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==1,]$match.ratio, 2), sep = "")
      neg.91.to.neg.52 <- rbind(neg.91.to.neg.52, paste(round(elec.data[elec.data$week==1,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==1,]$widower.effect.est / round(elec.data[elec.data$week==1,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      neg.91.to.neg.52 <- rbind(neg.91.to.neg.52, paste(" (", round(elec.data[elec.data$week==1,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==1,]$widower.effect.upper, 3), ")", sep = ""))
      neg.52.to.neg.15 <- paste("n = ", elec.data[elec.data$week==2,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==2,]$match.ratio, 2), sep = "")
      neg.52.to.neg.15 <- rbind(neg.52.to.neg.15, paste(round(elec.data[elec.data$week==2,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==2,]$widower.effect.est / round(elec.data[elec.data$week==2,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      neg.52.to.neg.15 <- rbind(neg.52.to.neg.15, paste(" (", round(elec.data[elec.data$week==2,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==2,]$widower.effect.upper, 3), ")", sep = ""))
      neg.15.to.0 <- paste("n = ", elec.data[elec.data$week==3,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==3,]$match.ratio, 2), sep = "")
      neg.15.to.0 <- rbind(neg.15.to.0, paste(round(elec.data[elec.data$week==3,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==3,]$widower.effect.est / round(elec.data[elec.data$week==3,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      neg.15.to.0 <- rbind(neg.15.to.0, paste(" (", round(elec.data[elec.data$week==3,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==3,]$widower.effect.upper, 3), ")", sep = ""))

      elec.vector <- cbind(neg.91.to.neg.52, neg.52.to.neg.15, neg.15.to.0, blank, blank, blank)
      assign(paste(option.name[i], election.name[m], sep = " - "), elec.vector[c(2, 3, 1), 1:6])

      m <- 2
      elec.data <- get(file[n + m])

      neg.15.to.0 <- paste("n = ", elec.data[elec.data$week==3,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==3,]$match.ratio, 2), sep = "")
      neg.15.to.0 <- rbind(neg.15.to.0, paste(round(elec.data[elec.data$week==3,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==3,]$widower.effect.est / round(elec.data[elec.data$week==3,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      neg.15.to.0 <- rbind(neg.15.to.0, paste(" (", round(elec.data[elec.data$week==3,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==3,]$widower.effect.upper, 3), ")", sep = ""))
      x.0.to.15 <- paste("n = ", elec.data[elec.data$week==5,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==5,]$match.ratio, 2), sep = "")
      x.0.to.15 <- rbind(x.0.to.15, paste(round(elec.data[elec.data$week==5,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==5,]$widower.effect.est / round(elec.data[elec.data$week==5,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      x.0.to.15 <- rbind(x.0.to.15, paste(" (", round(elec.data[elec.data$week==5,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==5,]$widower.effect.upper, 3), ")", sep = ""))
      x.15.to.52 <- paste("n = ", elec.data[elec.data$week==6,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==6,]$match.ratio, 2), sep = "")
      x.15.to.52 <- rbind(x.15.to.52, paste(round(elec.data[elec.data$week==6,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==6,]$widower.effect.est / round(elec.data[elec.data$week==6,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      x.15.to.52 <- rbind(x.15.to.52, paste(" (", round(elec.data[elec.data$week==6,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==6,]$widower.effect.upper, 3), ")", sep = ""))

      elec.vector <- cbind(blank, blank, neg.15.to.0, x.0.to.15, x.15.to.52, blank)
      assign(paste(option.name[i], election.name[m], sep = " - "), elec.vector[c(2, 3, 1), 1:6])

      m <- 3
      elec.data <- get(file[n + m])

      neg.15.to.0 <- paste("n = ", elec.data[elec.data$week==3,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==3,]$match.ratio, 2), sep = "")
      neg.15.to.0 <- rbind(neg.15.to.0, paste(round(elec.data[elec.data$week==3,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==3,]$widower.effect.est / round(elec.data[elec.data$week==3,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      neg.15.to.0 <- rbind(neg.15.to.0, paste(" (", round(elec.data[elec.data$week==3,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==3,]$widower.effect.upper, 3), ")", sep = ""))
      x.0.to.15 <- paste("n = ", elec.data[elec.data$week==5,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==5,]$match.ratio, 2), sep = "")
      x.0.to.15 <- rbind(x.0.to.15, paste(round(elec.data[elec.data$week==5,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==5,]$widower.effect.est / round(elec.data[elec.data$week==5,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      x.0.to.15 <- rbind(x.0.to.15, paste(" (", round(elec.data[elec.data$week==5,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==5,]$widower.effect.upper, 3), ")", sep = ""))
      x.15.to.52 <- paste("n = ", elec.data[elec.data$week==6,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==6,]$match.ratio, 2), sep = "")
      x.15.to.52 <- rbind(x.15.to.52, paste(round(elec.data[elec.data$week==6,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==6,]$widower.effect.est / round(elec.data[elec.data$week==6,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      x.15.to.52 <- rbind(x.15.to.52, paste(" (", round(elec.data[elec.data$week==6,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==6,]$widower.effect.upper, 3), ")", sep = ""))
      x.52.to.80 <- paste("n = ", elec.data[elec.data$week==7,]$cases.inweek, ", m = ", round(elec.data[elec.data$week==7,]$match.ratio, 2), sep = "")
      x.52.to.80 <- rbind(x.52.to.80, paste(round(elec.data[elec.data$week==7,]$widower.effect.est, 3)," (", 100 * round(elec.data[elec.data$week==7,]$widower.effect.est / round(elec.data[elec.data$week==7,]$control.turnout.est, 3), 3), "\\%)", sep = ""))
      x.52.to.80 <- rbind(x.52.to.80, paste(" (", round(elec.data[elec.data$week==7,]$widower.effect.lower, 3), ", ", round(elec.data[elec.data$week==7,]$widower.effect.upper, 3), ")", sep = ""))

      elec.vector <- cbind(blank, blank, neg.15.to.0, x.0.to.15, x.15.to.52, x.52.to.80)
      assign(paste(option.name[i], election.name[m], sep = " - "), elec.vector[c(2, 3, 1), 1:6])

      table <- rbind(table, cbind("", "", "", "", "", ""), get(paste(option.name[i], election.name[1], sep = " - ")), get(paste(option.name[i], election.name[2], sep = " - ")), get(paste(option.name[i], election.name[3], sep = " - ")))
      rownames(table)[(10*(i - 1) + 1):(10*(i - 1) + 10)] <- c(paste("*", option.name[i], "*", sep = ""), paste("    /", election.name[1], "/", sep = ""), "", "", paste("    /", election.name[2], "/", sep = ""), "", "", paste("    /", election.name[3], "/", sep = ""), "", "")
      }

  colnames(table) <- column.names

  write.table(table, "plots/Aggregate_table.csv", quote = FALSE, sep = "|")
